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ABSTRACT 



We investigate non-spherical behavior of gas accreting onto a central supermas- 
sive black hole. Assuming optically thin conditions, we include radiative cooling and 
radiative heating by the central X-ray source. Our simulations are performed using 
the 3D SPH code GADGET-3 and are compared to theoretical predictions as well as 
to ID simulations performed using the grid code ZEUS. As found in earlier ID studies, 
our 3D simulations show that the accretion mode depends on the X-ray luminosity 
{Lx) for a fixed density at infinity and accretion efficiency. In the low Lx limit, gas 
accretes in a stable, spherically symmetric fashion. In the high Lx limit, the inner gas 
is significantly heated up and expands, reducing the central mass inflow rate. The ex- 
panding gas can turn into a strong enough outflow capable of expelling most of the gas 
at larger radii. For some intermediate Lx: the accretion flow becomes unstable devel- 
oping prominent non-spherical features. Our detailed analysis and tests show that the 
key reason for this unstable non-spherical nature of the flow is thermal instability (TI). 
Small perturbations of the initially spherically symmetric accretion flow that is heated 
by the intermediate Lx quickly grow to form cold and dense clumps surrounded by 
overheated low density regions. The cold clumps continue their inward motion forming 
filamentary structures; while the hot infalling gas slows down because of buoyancy and 
can even start outflowing through the channels in between the filaments. We measured 
various local and global properties of our solutions. In particular, we found that the 
ratio between the mass inflow rates of the cold and hot gas is a dynamical quantity 
depending on several factors: time, spatial location, and Lx] and ranges between 
and 4. We briefly discuss astrophysical implications of such Tl-driven fragmentation 
of accreting gas on the formation of clouds in narrow and broad line regions of AGN, 
the formation of stars, and the observed variability of the AGN luminiosity. 

Key words: accretion — galaxies: nuclei — hydrodynamics — methods: numerical 
— radiation mechanisms: general. 



1 INTRODUCTION 

Active galaxies are believed to be powered by accretion of 
matte r onto th e central supcrmassivc black holes (SMBHs) 
(e.g.. ISalpeterl [1964: Lvndcn-BcU 1969; Bkndford 1979: 
Ozernoi fc Reinhardtl Il978: Bal ick fc HeckmanI Il98d: iReed 



1984 



Fabian fc Crawford ,1990r iKormendv fc Richstone 



iggalKauffmann fc Haehneltll2000l : fFerrarese fc Fordl2005h . 



Active Galactic Nuclei (AGN) are argued to influence cos- 
mological galaxy formation in the form of feedback, whereby 
the overall pr operties of a galaxy can be regulated by its cen- 
tral BH (e.g.. lSilk fc Reesill998l : [Kin3l200a : IWvithe fc Loebl 
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20031 : Granato et al. 20041: Murray. Qua.taert fc Thompson 



2005'; Begelman & Nath"2005': 'Bcst"2007VCiotti fc Ostriker 



i 2007c .Pipino. Silk fc Matteucci 2009 ; Ostriker et al.l I2OI0I : 
iDi Matteo et al.ll2011 ) . Observationally AGN are character- 
ized as a compact continuum source around a central SMBH, 
surrounded by a much larger emission line region consisting 
of multiphase gas clouds. The continuum emission shows a 
relatively flat spectral energy distribution over a broad wave- 
band from radio to gamma rays, and is modeled as coming 
from the central engine composed of a gas accretion disk of 
size 10"'' - 10"^ pc. The broad-line region (BLR) extends 
up to 0.01 — 1 pc and is composed of higher-density, higher- 
velocity gas; while the outer narrow-line regions (NLR) are 
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considered to be located at distances up to 10 — 1000 pc 
from the central SMBH. 

Modeling gas accretion onto SMBHs and resulting AGN 
feedback in a cosmological context is computationally chal- 
lenging, because a large dynamic range of length scales is 
involved: BH accretion on sub-pc scale, to galaxy physics 
on hundreds of kpc scale. Therefore SMBH accretion and 
feedback is usually t reated by a sub-grid model in cosmologi- 
cal si mulations fe.g. jDi Matteo et al.ll200i : lBooth fc Schavd 
I2OO9I ). At the same time, gas accretion flow within the Bondi 
radius of a SMBH has ju st started to be resolved in observa- 
tions (|Wong et alj|201ll l. Hence numerical studies resolving 
the Bondi radius are required for a meaningful comparison, 
and it would be useful to improve a sub-grid model of AGN 
feedback for cosmological simul ations based on t he results of 
small-scale simulations. In Bara i. Proga fc Naga minc (2 01ll . 
hereafter Paper I), we studied spherically-symmetric Bondi 
accretion of gas onto a SMBH, and started to explore the 
effects of radiative heating and cooling. 

A number of recent studies resolving the Bondi ra- 
dius have focused on accretion onto intermediate-mass 
BHs regulated by feedback from X-ray and UV radia- 
t ion emitted isotropicall y near the BH. Analytical work 

( Milosavlievic et al. |2009|) as well as ID and 2D simulations 
Milosavlicvi c. Couch fc Brommll2009l : iPark fc Ric"ottill201ll . 
[2012 ) find that photoionization heating and radiation pres- 
sure significantly reduce the steady-state accretion rate to 
a fraction of the Eddington rate, or over 2 orders of magni- 
tude b elow the Bondi rate. However, in ID simulations by,Li| 
l|201lll the self-gravity of the gas eventually overcomes the 
radiative feedback effects, and boosts the accretion to the 
Eddington r ate after one free-fall timescale. IPark fc Ricottil 
(|201ll . T201^ ') also find the development of Rayleigh- Taylor 
instabilities in their 2D simulations. In these work the BH 
luminosity is made to be proportional to the rate at which 
gas accretes across the inner boundary, whereas in our study 
the central luminosity has a constant value in each simula- 
tion because we want to have a simplified setup with a non- 
variable source to enable us cleanly investigate the physical 
processes influencing the accretion flow. 

In Paper I we discovered indications of thermal in- 
stability (TI) in some runs, which we investigate in de- 
tail in the current work. Studying TI in a dynamical flow 
in 3D is challenging, because the effects of a gravitational 
field significantly com plicate the development of TI (e.g., 
iBalbus fc Sokeilll989h . 

The foundations for understanding astrophysical ther- 
mal i nstability were laid in the classical paper by iFieldl 
The general criterion for TI to occur is: 



dc 

dS 



<0, 



(1) 



where C is the net cooling-heating rate, S is the specific en- 
tropy, and ^ is a thermodynamic variable that is kept con- 
stant for a given perturbation. Subsequently the theory of TI 
has been studied over decades, including analytical develop- 
ments: local TI in static and dynamical sy stems (e.g., JBalbus 
1 19861 : iBalbus fc So"ke3ll989l : lBalbuslll995l). global TI b y lin- 
ear perturbation analysis (e.g., Yamada fc Nishill200ll ). and 
numerical hydrodynamical calculations to investigate the 
pertu rbations in the non-linear regime (e.g.. lHattori fc Habd 
Il990l) . 



TI has been invoked to explain observed features in 
v arious astrophysical domains: solar prominen ces (e.g., 
ParkeJ Il953l : iKarpen. Picone fc Dahlburd Il988l l: forma- 
tion of cold gas clouds and cl umpy substructure in the 



inter ste llar medi um i^-Z-j Schwarz. McCrav fc SteirJ 



I972I: IParravand 119871 : iHennebelle fc Perauld 



Burkert fc LinI 2000l: Kovama fc Inutsukal 



Sanchez-Salcedo. Vazguez-Semadeni fc Gazol 



199E 



2002 



2002 



Inutsuka. Kovama fc Inouel I2OO5I ); origin of gl obular clus- 



ters during the collapse of a protogal axy (e.g., Fall fc ReesI 
19850; form at^ion of galaxies (e.g., iGold fc Hovld ~959l : 
Sofug [l969l ): general cooling flows in galax i es ari d galaxy 
clust er cores (e.g., Mathews fc Bregman 19781: iNulsenl 



19861: iMalagoh. Rosner fc Bodd 



Yoshida. Hattori fc Habd Il99ll : IGuo. Oh fc Ruszkowskl 
20081 ): the multi-phase structure of cold filaments con- 



1987 



MeiksinI 



1988 



densing out of the hot int racluster medium in galaxy 
clusters and group s (e.g., [Sh arma. Parrish fc Quataerg 
I2OIOI : iMcCourt et"aLll201ll ) . 

More relevant for our work are the applications of TI 
in black hole accretion systems. Spherically symmetric ac- 
cretion flow onto a compact object is preheated b y radi - 
ation from the central X-ray source. lOstriker et al.l (Il976l ) 
showed the preheating suppresses accretion, such that above 
a X-ray luminosity 0.00 5 of the Eddington lim it steady flow 
is impossible. iCowie, Ostriker fc Starlj (|l978l ) found time- 
dependent behavior of the preheated accretion for a much 
wider range of sourc e parameters, and the presence of in - 
stabilities. However, iBisnovatvi-Kogan fc Blinnikovl (Il980l ) 
found that numerical solutions exist for any luminosity be- 
low the Eddington limit, but they did not investigate the 
stab ility of the so l ution s in detail. Performing stability anal- 
yses ISteUingwerj (| 19821 ) found globally stable solutions, but 
those become subject to a drift instability at luminosities 
over the Eddingt on limit exce eding ~ 0.03. Hydrodynami- 
cal simulations bv lProgal (|2007l ) show that with an accretion 
luminosity of 0.6 of the Eddington value, the gas flow settles 
into a steady state with two components: an equatorial in- 
flow, and a bipolar inflow/outflow with the outflow leaving 
the system along the disk rotational axis. There the tem- 
perature at the outer radius was flxed at a relatively high 
value (2 x lO'^ K), so the gas did not go throug h a ther- 
mally unstable phase. Also lPark fc Ricottil (|2012 h flnd that 
a quasi-steady BH accretion is possible at the Eddington 
luminosity when the gas density is above a critical value 
inversely proportional to the BH mass. 

Investigating sphe rical accretion of gas irra diated by a 
quasar-like continuum, iKrolik fc LondonI (|l983*) found spe- 
ciflc criteria for when TI may disrupt single-phase steady 
flow, in the lu minosity - effic iency (of mass to energy con- 
version) plane. iKrolikI I 19881 ) argued that broad emission 
line clouds in AGN originates in cooler condensations form- 
ing within a hot medium via finite-amplitude TI. How- 
ever, iMathews fc Doand l|l990l ) concluded that broad line- 
emitting clouds are unlikely to arise from thermally unstable 
condensations in a hotter medium, but are likely composed 
of higher-density ga s ejected from nearby sta rs or from a 
dense accretion disk. IWang. Cheng fc Lil (|2012l ) investigated 
the dynamics of clumps, considered to exist by thermal in- 
stability, embedded in and confined by advection-dominated 
accretion flows onto black holes. 

We flnd that a spherically accreting gas distribution 
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acted upon by a central gravitational potential, undergo- 
ing radiative heating and cooling, can be thermally unsta- 
ble, cool non-spherically and become clumpy. Hence results 
from our simulations show that fragmentation can occur in 
a simple system with minimal physical processes. Our study 
can give insight into the formation and evolution of clouds 
near AGN. TI operating in such a way can be responsible for 
the origin of two phase medium in the NLR, support star- 
formation near AGN, contribute to the observed variability 
of the luminosity. 

This paper is organized as follows. We describe our sim- 
ulations in We present and explain the results in 331 and 
discuss in ^ We summarize and conclude in ^ 



2 SIMULATIONS 

In Paper I we performed simulations of gas accretion onto a 
SMBH of mass 10* Af© within the radial range of 0.1-200 pc, 
using the 3D S PH code GADGET-3 (originally described by 
ISpringel|[2005l ). The gas experiences gravitational potential 
of the central SMBH only, and there is no self-gravity. We 
first tested how accurately the GADGET code can follow 
the Bondi problem, aiming to know what was required for 
it to reproduce the Bondi solution. Our simulations could 
reproduce adiabatic Bondi accretion for durations of a few 
dynamical times at the Bondi radius, and for longer times 
if the outer radius is increased. Next we investigated one 
of the modes of AGN feedback: radiation heating, by per- 
forming simulations of gas accretion without and with radia- 
tive heating by a central X-ray corona and radiative cooling. 
We found that artificial viscosity causes excessive unphysi- 
cal heating near the inner radius. We also observed feedback 
due to radiative heating in our simulations, as increasing the 
X-ray luminosity {Lx ) produced a lower central mass inflow 
rate, and an outflow developed for sufficiently high Lx- 

All the results of gas properties presented in Paper I 
were spherically-symmetric. However, for certain parame- 
ters, we observe that the spherical-symmetry is broken, and 
the flow becomes inhomogeneous and multi-phase. Here, we 
analyze these non-spherical simulations and present corre- 
sponding results. We investigate the gas behavior in detail, 
aiming to comprehend the underlying physics versus numer- 
ics at play. 

The radiative heating and cooling prescriptions used in 
our simulations are described in §2.4 of Paper I. Radiation 
from the central SMBH is considered to be in the form of 
a spherical X-ray emitting corona of luminosity Lx which 
irradiates the accretion flow, and we consider an optically 
thin case. The local X-ray radiation flux at a distance r 



from the central source is Fx ~ 
parameter ^ is defined as 



The photoionization 



(2) 



Here n — p/(/imp) is the local number density of gas, com- 
puted using the mass density p, mean molecular weight /i, 
and the proto n mass rrip. The h eating-cooling fu nction is 
adopt ed from iProga et al. I l|2000l ). originally from iBlondinI 
l|l994h . These are approximate analytic formulae for the 
heating and cooling rates of an X-ray irradiated optically- 
thin, cosmic abundance (in this context meaning close to so- 
lar metallicity) gas illuminated by a 10 keV bremsstrahlung 



spectrum. Radiation pressure of the gas has not been in- 
cluded. The net heating-cooling rate, £, is given by 

pc = (Gcompton + Gx - Lt,i) = Function of (5, T) , (3) 

which includes Compton heating and cooling. X-ray 
photoionization heating and recombination cooling, and 
bremsstrahlung and metal line cooling. The characteristic 
temperature of the bremsstrahlung radiation has a constant 
value Tx ~ 1.16 x 10** K. The heating-cooling rate depends 
on density, since ^ is a function of p and r. 

We compute an additional quantity for the analysis pre- 
sente d in this work: the pressure- ionization parameter H 
(e.g., iKrolik. McKee fc Tartedll98ll ). which is defined as 



Fx 
cP 



4:Tvr^cP ' 



(4) 



where c is the speed of light, and P is the gas pressure. 

The functional form of radiative heating and cooling in 
Equa tion ([3l) has been implemented in the gr id code ZEUS 
(e.g., iProgaTbOOTl : [k urosawa fc Proga|[2009al lbl). and there 
are many resulting insights from ZEUS simulations about 
the dynamics of accreting gas irradiated by a quasar. Con- 
sidering radiation field from a UV-emitting standard accre- 
tion disk alo ng with a sph erical central X-ray source and 
line driving, iProgal l|2007l ) found that the accretion flow 
settles into a steady state and ha s two components ([JT]). 
jProga, Ostriker & Kurosawa ('2008') added gas rotation (also 
in Kuro sawa fc Proga 2009a b). and saw an increase of the 
outward thermal energy flux, leading to fragmentation and 
time variability of the outflow. 

As a new feature, in Paper I we increased the outer 
radius to rout ~ 200 pc. This in turn all ows us to follow the 
evolution of gas that is cooler than in [Kurosawa fc Progal 
(|2009bl ) (rout ~ 10 pc), so that the conditions are met for 
the gas to undergo thermal instability. The gas near rout 
is subsonic and lies in the cool, thermally stable branch of 
the T — ^ radiative equilibrium curve f i|3.1.1|) . As it accretes 
inward it goes to the thermally unstable branch. 

Table [T] gives a list of four simulations presented here. 
The run numbers correspond to those in Paper I. The only 
parameter varied between the runs is the X-ray luminos- 
ity Lx, in units of the Eddington luminosity LEdd- These 
simulations are done by restarting the Run 23 (presented in 
Paper I, with Lx = 5x 10~^) at t = 1.4 Myr, and performing 
subsequent Runs 25 - 28 with Lx = 0.005,0.01,0.02,0.05. 
Some of the common parameter values are: 7run ~ 5/3, 
rin = 0.1 pc, rout = 200 pc, poo = 10~^^ g/cm^, and 
Too ~ 10^ K. The temperature at the outer radius corre- 
sponds to a Bondi radius of Rb = 183.9 pc. The Bondi 
accretion rate is Ms(7 = 5/3, poo, Too) = 4.9 x 10^^ g/s. 
This is 0.35 times the Eddington accretion rate, assuming a 
radiative efflciency of 0.1. 

An assumption of our work is a constant value of Lx in 
a simulation. This enables us to have a simple setup in order 
to investigate the physical origin of the non-spherical behav- 
ior of the accretion flow. A more realistic scenario would be 
Lx as a function of the mass accretion rate at the inner 
boundary. However this would lead to an additional depen- 
dence, which we leave for future studies. 
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3 RESULTS 

3.1 Lx = 0.01: Cold Clumps 

There is fragmentation and clumping in Run 26, as the gas 
cools and becomes denser in certain regions. The free-fall 
time at the outer radius is 2.98 Myr. Non-spherical features 
start to develop from ~ 1.6 Myr continuing as long as the 
simulation is run (2.12 Myr), occurring at r < 30 pc. 

Figure[T]shows the inner 30 pc cross-section in the [x— y] 
plane through z = aX t = 2.047 Myr, with the gas den- 
sity in the top-left, temperature in the top-right, photoion- 
ization parameter in the bottom-left, and Mach number in 
the bottom-right panels, overplotted with velocity vector ar- 
rows (the scale is indicated at the top- left in units of km/s). 
The figure illustrates the filamentary structure of the cold 
infiow. The cold-dense filaments at r > 5 pc have a lower- 
photoionization, and extend up to r ~ 20 — 30 pc. Hotter, 
lighter gas exists in-between the filaments. Both the hot and 
cold components fiow inward initially. The hot phase is sub- 
sonic, while the cold phase is supersonic reaching Mach ~ 10 
all the way down to rin. At t > 1.8 Myr some of the hot- 
ter gas between r ~ 25 — 30 pc start to flow outward, but 
there is no well-developed outflow, since the hot gas encoun- 
ters colder inflowing gas at larger radii. The photoionization 
parameter is high near the center despite the high central 
density. 

The details of the inner-most cold clumps are plotted in 
Figure [2l It is a cross-section of the inner 4 pc of the [x — y\ 
plane aX t = 2.047 Myr (same plane as Figure [T] zoomed- in 
here). However note that this cross-section uses a different 
color scheme and it has been plotted without the velocity 
vectors, in order to show the small-scale features clearly. 
The central r < 1 pc region is denser than the surrounding, 
since the stable free-fall gas density increases going inward 
(e.g.. Figure 7 in Paper I). In terms of other properties this 
central region is characterized by high temperature, high 
photoionization parameter, and supersonic velocity. 

We investigate the motion of the outer clumps as 
they fall into rin by making movie^ of multiple snapshots 
within 30, 12 and 4 pc cross-sections. After formation from 
thermally- unstable gas (details in tj3.1.1|) . as the clumps fall 
in they undergo a dynamic evolution. Some stretch and 
fragment producing two or more sub-clumps, caused by the 
gravitational held of the central SMBH. Sometimes smaller 
clumps merge with each other creating a bigger one. 

In the inner volume, at r < 0.5 pc, although the den- 
sity of the falling clumps remains significantly higher than 
the mean density, the temperature becomes comparable to 
the mean value at a given radius. Consequently we can see 
a dense clump in the density-plane, but it is not cold any 
more in the temperature-plane. An example in Figure [5] is 
the distinct dense blob at j/ ~ —0.5 pc, which is also visible 
as a low-photoionization and high-Mach clump. The cen- 
tral heating at r < 1 pc is mostly by adiabatic compression 
where the entropy of a gas particle remains constant. Radia- 
tive heating would dominate at a few pc for some particles. 
We have checked the role of artificial viscosity (given our 
results in Paper I), and found that although it increases 
in magnitude toward rin, it always remains sub-dominant 

^ http: / / adIibitum.oats.inaf.it /barai /AUPages / visualization.html 



compared to adiabatic heating in the simulations presented 
here. 

The gas particle properties as a function of radius at 
t — 2.0 Myr is plotted in Figure O It can be compared to 
Figure 15 in Paper 1 which shows the same simulation run 
at an earlier time, t — 1.5 Myr, where the flow is still stable 
and spherical. 

Figure O shows radial scatters at r < 30 pc as large 
as 2—3 orders of magnitude. It corresponds to the multi- 
phase structure of the gas, with a group of particles in the 
hot branch (T > 10^ K), another group in the cold branch 
(T < 3 X 10'' K), and some in transition between the two. 
The cold particles are adiabatically heated up at r < 0.7 — 1 
pc, with no particles below 10^ K at those inner radii. The 
temperature increases to > 10^ K near nn. Some particles 
are in free-fall, some are slower forming a scatter by a factor 
of a few, and some outflowing {vr > 0) between 7 — 25 pc. In 
the temperature versus ionization parameter planes {T — ^ 
and r — H), the particles are scattered around the respec- 
tive unstable branch (see tj3.1.ip of the radiative equilibrium 
curve. The group of particles departing from the radiative 
equilibrium curve with little scatter in the lower-right por- 
tion of these planes corresponds to the outflow near rout. 



3.1.1 Thermal Instability 

The slope of a perturbation in the [T — ^] and [T — H] planes 
can be used to determine if the gas is stable against such a 
perturbation. The relations for three common processes are: 

r oc ^^T (X H"^(adiabatic free fall, or, in Eq. Q A = S), 
T oc 5, H = constant (constant pressure, or, A — p), 
^ = constant, T oc (constant density, or, A = p). 

These relationships are shown as the blue (dotted, dashed, 
and dash-dotted) lines in the bottom row of Figure [l] 

According to our radiative heating and cooling prescrip- 
tions, the upper-left half of the Trad — £, equilibrium curve 
(red line in the bottom-left panel of Figures [3] and [3]) is the 
net cooling zone, and the lower-right half is the net heating 
zone. Let us consider a perturbation with constant pressure 
on a point on the T^ad — £, equilibrium curve. The accretion 
flow is stable in the parameter region where the instanta- 
neous slope of the equilibrium curve is shallower than T cx ^, 
and vice versa. The T^ad — £, equilibrium curve has the steep- 
est slope between ^ ~ 100 — 500, therefore it is easier to flnd 
a perturbation with a shallower slope than the equilibrium 
curve, making it most prone to instability there. We call 
this region as the unstable region. In our simulations we in- 
deed find that gas particles undergo TI within the unstable 
^-range, as described below. 

We also compute T — ^ solution of th e isobar i c insta bil- 
ity criterion given in the Equation (5) of lBalbr3 l|l986h for 
local TI in dynamical systems. Using ^ oc 1/p, the relevant 
expression becomes 




-1 = 0, 



which gives the green line in the bottom-left panel of Fig- 
ure |4l It exactly overlaps with the Trad — ? radiative equi- 
librium curve in the unstable region. 
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By selecting a single particle we study the time evolu- 
tion of its various radial properties, and its thermal history 
in the [T — ^] and [T — H] planes. Using this approach for 
particles selected from different parts of the [T — ^] plane 
(bottom- left panel in Figure |3]), we find that at initial times 
{t < 1.6 Myr, when the flow is still stable at all-r) every par- 
ticle has an increase in temperature while moving inward. 
After 1.6 Myr we find two broad categories of evolution. 
In the first kind, a particle cools down to T ~ 3 x 10* K, 
might undergo alternative periods of heating and cooling 
while continuing to move in, and finally heats up near rin. 
These cooled particles form the dense filaments seen in Fig- 
ures [TJ [2] In the second kind, a particle would never cool, 
but continue to heat up all the way to the center. These hot 
particles form the lighter gas in-between dense filaments. 

We argue that this cooling (which occurs in a non 
spherically-symmetric manner) is caused by the TI in accret- 
ing gas acted upon by the adiabatic and radiative processes. 
Investigating the thermal history, we find that a single parti- 
cle goes through both adiabatic and radiative process, each 
of which dominating at different times. We hereafter call 
these epochs as adiabatic and radiative regimes. In the ra- 
diative regime the particle follows the radiative equilibrium 
temperature (Trad), and in the adiabatic regime it is heated 
up following the adiabatic law (Tfl.a oc r~^). We determine 
the dominant process by comparing the terms in the du/dt 
(time derivative of the specific internal energy) equation, 
as well as comparing the particle evolution to the free-fall 
slopes in the [T — r] and [p — r] planes. We also verify that 
its entropy remains constant when adiabatic processes dom- 
inate. 

We note that the adiabatic free-fall relation, Tff,a oc 
has the same slope as the radiative equilibrium curve 
in the unstable region. Also in the [T— H] plane, the adiabatic 
free-fall relation (Tff ,a c>c Hg'^) has the same slope as Trad in 
the unstable region, where dT/dS. < (isobaric case). This 
can be seen in the bottom two panels of Figure H] as the 
overlapping red and blue-dotted lines. 

During the linear phase of TI growth, a subsonic par- 
ticle lying on an unstable part of the Trad — C equilibrium 
curve would move along the constant-pressure (T oc ^) slope. 
However our system is very dynamical and the non-linear de- 
velopment of TI can be different from classical descriptions. 
Therefore it might not end up on another stable equilibrium 
branch. The adiabatic heating (or, cooling) arising from in- 
ward (or, outward) gas motion can prevent a particle from 
following a path with constant-pressure all the way to the 
next stable position. 

We detect a distinct mode of cooling in the [T — ^] 
plane, which occurs as below. A particle while being in the 
adiabatic regime moves up to higher-T along the unstable 
part of Trad — This happens because of the equality of the 
adiabatic and radiative slopes in the unstable regions, as de- 
scribed before. At a later time the particle has a transition 
from adiabatic to radiative regime. Then TI operates, and 
the particle undergoes catastrophic cooling or heating (de- 
pending on which is dominant), and moves away from the 
unstable Trad — C part. It finally settles on one of the stable 
branches where T^ad — C has a smaller slope. 

Figure U] plots the evolution of a typical particle which 
goes through such a cooling instability, and demonstrate the 
above arguments. This particle undergoes evolution of the 



first kind: heats up initially (A— s-B-^C), cools (C to D), 
and finally heats up again near the center (D— >-E— >-F— >-G). 
In between there could be several episodes of alternative 
heating and cooling each lasting for short durations, as we 
describe below. 

Let us follow the particle evolution on the T^ad — £, plane 
in more detail. Starting from the IC at t = 1.4 Myr when 
it is at point A, r = 53 pc indicated by the triangle in each 
panel, the symbols denote uniform time intervals of 

0.004 Myr (except the two inner-most points F and G which 
are separated by ^ 0.001 Myr). The ending point G is at 
t = 1.8 Myr and r = 0.99 pc denoted by the square symbol 
in each panel, after which the particle is accreted within 
< 0.001 Myr. The IC has been taken from a run with a 
lower-Lx , but with the same density, which corresponded to 
a lower-^. On starting this run with a higher-ix the particle 
has a higher-^. Therefore it is radiatively heated up to Trad 
initially from point A to B, a phenomenon that happens for 
all the particles in the simulation. 

After this initial quick heating, the particle reaches the 
point B belonging to the unstable branch of T^ad- However 
by that time it enters the adiabatic regime at r ~ 50 pc, and 
moves up along the unstable Trad — C branch from point B 
to C due to adiabatic heating. This can be verified in the 
[T — r] and [p — r] planes where it follows the free-fall slope, 
and its entropy remains constant. At point C, r ~ 35 pc, it 
transits to the radiative regime with dominance of cooling 
and TI operates. It cools down rapidly from C to D deviating 
slightly to the left of equilibrium curve (deviation is more 
evident in the T — H plot), i.e. in the cooling zone. By point 
D, r ~ 15 pc, it reaches the cold-stable Trad - Then it becomes 
dominated by radiative heating for some time, and at r ~ 10 
pc undergoes another transition to adiabatic regime, all of 
which happens between points D and E. It reverts back and 
radiatively cools between r = 7 — 5 pc reaching point E. 
After that it heats up again while falling inward, first from 
E to F by radiative heating. The final heating from F to G 
at r < 2 pc is adiabatic, caused by the central compression 
of the gas. The transition of E— >-F^G can also be regarded 
as a TI, since it rapidly deviates from the stable branch of 
the equilibrium curve and landing on another stable branch. 



3.1.2 Timescale Comparison 

An important factor for TI to develop is the gas time scales, 
as ha s been discussed by other autho r s (e.g., Beltrametti 
198l| : iDavid. Bregman fc Se;ibl 1 19881 : iMathews fc Doane 
1990f ). If unstable perturbations occur in supersonic gas. 



then there are no resulting dynamical ch anges, since pres- 
sure forces play only a minor role (e.g., iKrolik fc LondonI 
119831) . In other words, the gas with thermal properties in 
the TI zone should be subsonic, otherwise there will not be 
enough time for TI to grow. For example this is the reason 
that simulations with Lx < 0.01 are smooth and spherically- 
symmetric. 

We investigate the origin of gas clumping in this run 
by comparing the relevant timescales. We compute the TI 
growth tim e (tgr) using the growth rate (ni) from Equa- 
tion 



(31) of iFieldl l|l965f ). which in our notation is 



1 

ni 



jkTo 



po 



dC 

dp 



(8) 
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Here po,To are the unperturbed density and temperature; 
where we use the current gas properties at a given radius 
when the flow is still spherically-symmetric. We analyze the 

ratio of tg,. with the free-fall timescale, t// = ( 2GMbh I ' 
which is plotted in Figure [5] top panel at three simulation 
epochs. We flnd positive ni in some radial ranges represent- 
ing regions where TI can grow, and negative ni elsewhere 
indicating regions stable to TI seen as the radial gaps in 
Figure \5\ 

We flnd that within the positive ni regions, tgr starts 
to become less than tff inside a range r ~ 10 — 80 pc soon 
after start of the simulation at t = 1.4 Myr. But the gas do 
not clump as soon and over that entire r-range because it is 
still supersonic, seen in the bottom panel of Figure [5] where 
Mach number is plotted. As discussed at the beginning of 
tj3.1l non-spherical features flrst appear from ~ 1.6 Myr, at 
10 — 20 pc. This is when and where the gas becomes subsonic 
is addition to tgr <tff, fro m Figure [5l and therefor e clump- 
ing develops, as expected. iKrolik fc LondonI (|l983l ) showed 
that unstable accretion flow, which is subsonic, can be dis- 
rupted when tgr /tff becomes small, exactly as we see in our 
simulations. 

As a quantitative analysis of the growth of instability, 
we plot the range of perturbation amplitude in Figure [G] It 
shows the time evolution of the temperature and density per- 
turbation at r = 10 pc within a width of Ar = 1 pc. Track- 
ing the gas particles lying within that radial bin, we plot 
the minimum, average, and maximum values as a function 
of time. The corresponding average value ai t — 1.42 Myr 
is considered as the unperturbed gas property. This time is 
close to the initial time of 1.4 Myr but distant enough such 
that the gas has adjusted to the effect of higher-Lx at the 
beginning, and when the flow is still stable and spherical. We 
plot the statistics relative to the initial, unperturbed value 
at t = 1.42 Myr. 

The minimum, maximum and average values are all 
close to the unperturbed average at t = 1.42 — 1.55 Myr, and 
they start to deviate after that. The perturbations show an 
exponential growth at first, and then come to a saturation 
for the remaining of the simulation. There is an exponen- 
tial decrease of Tmin between 1.55 — 1.65 Myr reaching a 
value 0.02, indicating the growth of cooling instability, and 
rmin remains almost constant after that. A gradual increase 
is seen for Tmax from 1.55 Myr to 2.17 Myr. The value of 
(T) decreases to 0.6 at t — 2.17 Myr. The value of pmax 
exponentially increases to 200 between t — 1.55 — 1.65 Myr, 
indicating the growth of cooled, dense clumps due to TI. It 
subsequently fluctuates between 100 — 200. Following a sim- 
ilar trend, (p) exponentially increases to 10 by f = 1.7 Myr, 
and remains constant after that. The value of pmin gradually 
decreases to 0.4 at t — 2.17 Myr. 

3.1.3 Fraction of Cold vs. Hot Gas in Accretion Flow 

In our simulations we have cold and hot phases of gas form- 
ing self-consistently in the accretion flow acted upon by ra- 
diative heating and cooling (whose functional form are in 
Paper I, Section 2.4) and undergoing TI. Furthermore we 
resolve the Bondi radius. The multi-phase gas has a temper- 
ature varying from 2 x 10* K to 10^ — 10^ K at the same r in 
Run 26. We choose a limiting temperature of Tum = 10^ K 



to separate the cold and hot phases. We measure the mass 
inflow rate of cold and hot gas (Min,coid and Min,hot) by 
counting the particles with Vr < 0, below and above Tum- 
We define the ratio of mass inflow rates of cold gas over hot: 
X = Min.coid/Min.hot- It is similar to the a-parameter multi- 
plied to the Bondi accretion rate in some of the galaxy sim- 
ulations (discussed in §[4]), which accounts for the fraction 
of cold, dense gas that is missed due to lack of resolution. 

We find that x is a- dynamical quantity depending on 
several factors: time, spatial location, and X-ray luminosity. 
Figure [7| shows the evolution of x in Run 26 as a func- 
tion of radius and time. The spatial distribution (in the left 
panel) indicates that x is non-zero only within certain r- 
range, where there is cold gas. In this run with Lx = 0.01, 
the cold gas starts appearing within r ~ 5 — 20 pc, where 
cooling instability occurs. The range increases to r ~ 1 — 30 
pc at later times. The Bondi radius of the hot gas (with 
T = 10® K) is ~ 18 pc. There is no cold gas within r < 1 pc 
due to the strong adiabatic heating near the center. The 
maximum value of x varies between 1 — 5. 

We also examined the temporal evolution of x 
r = lOpc by counting all particles whose smoothing kernel 
touches r = 10 pc. The result is plotted in the right panel of 
Figure [3 showing that x increases with time. The value of 
X increases quickly during the flrst 0.2 Myr, and after that 
it remains almost constant at a value x ~ 3 — 4. 

3.2 Lx = 0.02: Cold Clumps and Hot Bubbles 

Run 27 shows more complex non-spherical behavior in the 
form of cold clumps and buoyantly rising hot bubbles. 
The fragmentation and clumping start to develop from 
t ~ 1.7 Myr, and continues until the simulation ends at 
t = 2.46 Myr. The gas cools and becomes denser over 
r ~ 40 — 100 pc, following the physics of TI described in 
tj3.1.1l This radial range where Tl-induced fragmentation 
occurs is at larger radii than in Run 26 {Lx ~ 0.01). Be- 
cause with twice higher Lx in Run 27, the ^-range for TI 
(steepest part of the T^ad — (. equilibrium curve) is satisfled 
at larger radii. 

Figure [8] shows the time evolution of inner 100 pc cross- 
section of Run 27 in the [x — z] plane through y = 0. The 
gas is heated up and expands into a spherically-symmetric 
outflow in the inner volume (r < 40 pc), as the top row 
at t = 1.86 Myr illustrates. This initial outflow decreases 
the mass inflow rate at rin by a factor of ~ 1000 at t = 
1.4 — 2.2 Myr (discussed in i ]3.4|) . The surrounding gas at 
r = 40 — 70pc is inflowing, and has fragmented in certain 
regions. The central outflow collides with this clumpy inflow 
at r ~ 40 pc, creating a shock. It is visible as a ring-like 
density enhancement, which corresponds to the pile-up of 
gas outflown from the center and the cold clumps from outer 
regions trying to move in. 

The dense clumps forming at r ^ 40 pc continue to 
move inward, disrupting the spherical-symmetry of the cen- 
tral hotter outflow. By t — 2.12 Myr (middle row of Fig- 
ure [S]), the whole volume of r < 80 pc shows non-spherical 
features. We see the formation of a hot, less-dense bubble 
moving outward along the positive x-axis just off-center. 

The clumps finally start accreting into rin at t J5 
2.2 Myr. This causes the central mass infiow rate in this run 
to rise again. The bottom row (t = 2.46 Myr) of Figure [8] 
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shows an inhomogeneous multi-phase gas motion. There are 
signatures of a few hot bubbles buoyantly rising, which are 
visible as the density depressions in the positive-s plane. 
The bubbles are not well developed here as compared to 
Run 28 ( i)3.3|l . because during their propagation they en- 
counter clumpy medium and become distorted. 

3.3 Lx = 0.05: Hot Bubble 

The gas is heated up strongly by the highest Lx in Run 
28, and there is a spherical outflow. There is as well as an 
anisotropic hot buoyant bubble developing because of con- 
vective instability (timescales compared in H3.3.f p . A large 
fraction of gas is ejected out of the computational volume, 
and the central mass inflow rate drops drastically ( i]3.4|) . 
Convective instability becomes important in accretion flows 
locally dominate d by external heating, primarily in regions 
of subsonic flow jSalbus fc Sokeij|l989l ). Therefore this dy- 
namical instability is prominent in our larger Lx runs, where 
the gas is heated up and subsonic. 

A spherically-symmetric outflow develops inside of r ~ 
20 pc from t — f.46Myr, when the outer volume still has a 
spherical inflow. At i = 1.75 Myr there is still a 30 pc wide 
shell of inflowing gas at r ~ 140 pc. The gas is outflowing 
all over the volume after t — 1.8 Myr. 

Figure [9] shows the cross-section of the whole compu- 
tational volume ±200 pc of Run 28 in the [y — z] plane 
through a; = at two time epochs. There is the formation of 
a well-defined hot, less-dense bubble, which buoyantly rises 
from the center along the negative 2-axis. In Figure [5] it 
corresponds to the anisotropic jet-like structure: redder in 
temperature-plane, and bluer in density-plane. It starts at 
t ~ 1.7 Myr as a temperature enhancement and density de- 
pression at z ~ — Ipc. The top row a,t t = 1.8 Myr shows 
the bubble near its formation. It grows with time, with gas 
moving faster inside it than the surroundings. The bottom 
row at t = 3.0 Myr shows a developed stage of the bubble. 
It has a jet-like structure, being fed by hot gas near rin. 
Its head presents a bow shock-like morphology on top of a 
mushroom cloud; the interior gas collides with the relatively 
cold, ambient medium and produces a shock larger than the 
jet. The bubble head reaches rout at i ~ 3.8 Myr, and the 
rest of it dissipates after that. 

In this run the remaining of the gas outflow is 
spherically-symmetric. Propagating through relatively uni- 
form ambient gas, the bubble is hence well developed here, 
as compared to the distorted ones in Run 27. 

The properties of gas from two characteristic regions 
of the computational volume in this run at t = 2.0 Myr is 
plotted in Figure [TOl The cyan plus symbols are particles at 
r < 60 pc overlapping with the negative-y axis through their 
smoothing lengths, hence correspond to part of the bubble. 
The black plus symbols are particles overlapping with the 
positive-^ axis, representing the gas in the remaining vol- 
ume, and show spherically-symmetric behavior. 

The black points can be compared to Figure 16 in Pa- 
per I, which shows the same simulation run at an earlier time 
{t — 1.5 Myr) where the flow is still spherically-symmetric, 
and the spherical outflow had already developed. The inner 
gas in Figure 16 is pushed outward continuously, and the 
points are absent at r < 1 pc here in Figure 1101 

Outside the bubble the radial velocity shows mostly out- 



flowing gas, with some particles inflowing at r < 10 pc. The 
strong outflow in this run produces a density profile declin- 
ing with decreasing radius, which is an opposite radial slope 
compared to all other runs (e.g.. Figure [3] top-right panel). 
The temperature and entropy plots show increasing heating 
at smaller r, as expected by the high Lx- In the tempera- 
ture versus ionization parameter planes (r — ^ and T — H), 
the particles are near the hot, stable branch of the radia- 
tive equilibrium curve. But they are below T^ad since the 
outflowing gas is adiabatically cooling. 

The bubble (cyan points) has greater outward veloc- 
ity, lower density, higher temperature and higher entropy 
compared to the gas at the same radius. The gas inside the 
bubble has roughly constant density. It also has a roughly 
constant entropy, indicating that indeed the bubble is buoy- 
antly rising and convectively unstable. There is a slight de- 
crease in temperature with increasing radius, which is an 
indication of adiabatic cooling. 

The time evolution of a single particle inside the bub- 
ble is shown as the blue curve overplotted in each panel of 
Figure 1101 accompanied by the letters A— G in blue. The 
starting point A is indicated by the triangle, which is from 
the IC at t = 1.4 Myr and r = 11 pc. The ending point 
G denoted by the square is at t = 3.8 Myr when it is at 
r = 163 pc, and the direction of evolution is shown by the 
arrow. The particle is heated up from A to E at varying rates 
before reaching the bubble at point E, but finally cools in- 
side it from E to G. Apart from an initial rise from A to C, 
the density decreases monotonically from C to G at differ- 
ent rates. The entropy increases at first from A to E, and 
becomes almost constant inside the bubble from E to G. 

There is an initial quick radiative heating from A to 
B, since the IC has been taken from a run with a lower ^, 
after which the particle reaches an unstable branch of Trad 
at point B. It then undergoes heating instability from point 
B, because of the high-Lx in this run, and tries to reach 
the upper hot-stable branch. Once it becomes a part of the 
bubble at point E, its temperature remains almost constant 
as it is pushed out fast from E to F. It cools adiabatically 
while moving outward from F to G. The position of the 
particle at t = 2.0 Myr is denoted by the orange filled circle. 



3.3.1 Timescale Comparison 

The outflowing gas is always subsonic in this run, hence is 
prone to TI if the timescales are appropriate. We estimate 
the outflow motion timescale as: tflow = r/|-i;r|, and com- 
pare it with the TI growth time {tgr, Eq. |5}. The flow is 
stable initially (we find negative growth rate ni) at sim- 
ulation times ^ 1.6 Myr, and hence the gas do not show 
prominent clumping like in the previous runs. The flow be- 
comes unstable later, as we find tgr < tflow within a range 
r ^ 100— 120 pc occurring between a time interval ~ 1.6 — 2 
Myr. We argue that TI comes to play at that epoch lead- 
ing to the formation of a cold, dense spherical shell (which 
is clumpy at the beginning), which forms at ~ 100 pc and 
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moves outward. The shell can be seen in the top right panel 
of Figure [O] and clearly in the movi^f]. 

We verify the convective origin of the bubble by compar- 
ing the relevant timescales. We compute the Brunt- Vaisala 
timescale (tsv) for convective instability using the effec- 
tive Brunt- Vaisala freq uency (ujbv) from Equation (4.4) of 
iBalbus fc Sokeil (|l989D . 



UJBV 



p or J \ 



3 dlnP _ dlnp 
5 



dr 



dr 



-1/2 



(9) 



It is estimated when the flow is nearly steady and spherical, 
using the median values of gas density and pressure at each 
r. The ratio of tsv with the free-fall time is plotted in Fig- 
ure mi We find ujIv > in some radial ranges indicating 
convectively unstable gas. At the other radii we see convec- 
tively stable gas, where uJbv < Oi indicated by the radial 
gaps in Figure [TT] 

Within the unstable zone, we see that tsv becomes 
smaller than t// at some r which can be attributed to the 
development of the outflow and bubble in this run. In Fig- 
ure [TT] we find a convective instability point at r ~ 15 pc 
at simulation time 1.46 Myr where tsv < iff, which prop- 
agates outward with time. This represents the spherically 
symmetric outflow. There are also certain radii in the range 
(0.1 — 2) pc where tsv < tff, in between simulation times 
1.5 — 1.7 Myr. One such occurrence of convective instability 
is responsible for the development of the anisotropic bubble. 



3.4 Comparing All Runs 

We observe various modes of gas accretion for different val- 
ues of Lx- The global features of the flow in each run are 
summarized in Table [T] In Run 25 with Lx ~ 0.005, there 
is a stable inflow maintaining spherical-symmetry. This run 
is listed in the table for comparison, and has been dis- 
cussed in Paper I. The gas heats up as Lx is increased, 
and the accretion flow undergoes unstable motion causing 
non-spherical features. With Lx ~ 0.01 and 0.02 the gas 
cools and fragments forming multiple cold, dense clumps 
and attains a multi-phase structure. In the runs where a 
spherically-symmetric outflow is produced in the inner vol- 
ume {Lx = 0.02, 0.05) initially, there is also the formation 
of hot bubbles rising buoyantly from the center. 

Figure [12] shows the time evolution of the mass inflow 
rate at the inner radius (Min,ri„) for the four runs. We see 
that Min,ri„ drops at a varying rate as Lx is increased, be- 
cause the gas is heated up and expands. The inflow continues 
in the inner volume (r < 40 pc) with Lx ~ 0.005 and 0.01 at 
a rate 2 — 4 Mq/jt, with Min,ri„ reduced by a factor of 0.84 
in the latter run. When Lx ~ 0.02 or higher, a spherical 
outflow develops at r < 40 pc, because the gas is efficiently 
heated up by the central radiation. This outflow suppresses 
the accretion rate by 3 — 4 orders of magnitude between 
1.4 — 1.6 Myr in Runs 27 and 28. The mass inflow rate in 
Run 27 rises again after 2.2 Myr, because of the accretion 
of cooler clumps (details in i]3.2|) . There is however a strong 
outflow with Lx ~ 0.05 which expels a substantial fraction 



http:/ /adlibitum.oats.inaf.it/barai/AllPages/Images- 
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of gas out of the computational volume, and Min,ri„ never 
rises within 5 Myr. 

With Lx = 0.01 (Run 26) the hotter outflow reaches a 
radial velocity Vout ^ 300km/s between 10 — 30 pc. In Run 
27 {Lx ~ 0.02), initially there is a spherically-symmetric 
outflow with Vout ^ 200km/s at r < 40 pc. The outflow be- 
comes stronger at t ^ 2.0 Myr when hotter gas and buoyant 
bubbles become channelled in-between the inflowing cool, 
dense clumps, reaching Vout ^ 800 km/s. With Lx ~ 0.05 
in Run 28, the spherical outflow has Vout ^ 200 km/s, and 
the anisotropic hot buoyant bubble consist of faster mov- 
ing gas with Wout ^ 600 km/s. In runs with Lx = 0.02 
and 0.05 the fastest outflowing gas is subsonic at the be- 
ginning, and becomes supersonic later with Mach number 
of 2 — 3. At t ^ 1.8 Myr there is always some gas at 
Vout ~ 100 — 200 km/s which is supersonic. Initially most 
of the gas has Vout smaller than the escape velocity (wosc). 
At t > 2 Myr some gas between 40 — 120 pc flows out with 

V ^ Vcsc ■ 

The mass outflow rate (Mout) as a function of radius at 
a few specifled times, and its time evolution at a given radius 
for the three runs are plotted in Figure [TH] The outflow rate 
is low with Lx = 0.01, only reaching Mout ^ O.lMo/yr at 
r ~ 20 pc. It increases in the higher Lx runs, because of the 
development of strong outflows. With Lx ~ 0.02 it rises up 
to Mout ~ IMq/jt between r = 20 - 80 pc. For Lx = 0.05 
it reaches Mout = 5-10 Mq /yr at r = 40 - 100 pc. The 
resuhs in our GADGET runs with Lx = 0.01 and 0.02 are 
roughly in the same range as that obtained in the ZEUS 
counterpart sim ulations, which gave Mout ~ 0.1 — 1 Mq/yt 
l|Kurosawa fc Proga,2009ai '). 



4 DISCUSSIONS 

Our analysis of the hot versus cold gas fractions in the 
accretion flow ( i]3.1.3p can be used to constrain the a- 
parameter in the expression for SMBH accretion rate used 
in the sub-grid models of cosmological simulations. In 
such models, the Bondi-Hoyle-Lyttleton mass accretion rate 
(jHovle fc Lvttletonll 19391 : iBondi fc Hovi3ll944l : lBondilll952l ) 
inferred from simulations is multiplie d by the factor with 



a value oi the order a = 100 (e.g.. ILii Matteo et al.l IzUUal : 
[Dubois et aIll201Gl : ILusso fc Ciottill201Gl '). One of the reasons 
given for such high values of a is the artificially low densities 
and high temperatures of the ISM obtained by smoothing 
on the s cale of computational resolution at the location of 
the BH ([Booth fc S chavc 20091: [Johansson. Naab fc BurkertI 
[20091 : [Khalatvan et al...200a ). Another reason given is that 
a volume- average of the Bondi-rates for the (spatially unre- 
solved) cold and hot phases of the ISM recovers a value of a 
close to 100 ([Siiacki. Springel fc Haehne h 2009). Effectively 
the galaxy simulations need to assume that a large fraction 
of unresolved cold gas accretes onto the central SMBH, in 
order to grow the BHs to currently observed masses. 

However we do not find a factor as large as 100 for the 
cold gas fraction, as our x ^ 1 — 5. In our simulations we 
do not see a large amount of gas hidden in the cold phase. 
Other physical processes, like self-gravity and rotation which 
we do not have in our current simulations, might contribute 
to the conversion of hot accreting gas to the cold phase and 
increase x- The value of a should be a resolution dependent 
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quantity. A self-consistent determination of a needs multi- 
leveled high-resolution approach to study both star forma- 
tion and gas ac cretion onto SMBH (see recent attempt by 
iKim et al-feoill ). 

The mass inflow rate (Figure 1121 ^3.4p shows large 
fluctuations in Run 26 {Lx = 0.01) after 1.6 Myr, and in 
Run 27 (Lx = 0.02) after 2.3 Myr. This is because of the 
accretion of multiphase gas, with the spikes corresponding 
to the infall of co l d den se clumps. It is similar to the result 
by ISharma et al] (|201ll ). who found large variations of the 
mass accretion rate as a function of time when the cold TI- 
induced fllaments cross the inner boundary. The clumps in 



our simulations have density between 10 



10' 



j/cm , 



which overlap with the upper limit where NLR clouds are 
observed having electron density between lO'^ — lO'' cm~'^ 
(mass density 10~^^ — 10~^^ g/cva?). The temperature of the 
simulation clumps is between 10* — 10^ K at r ^ 1 pc, whose 
lower limit coincides with the observed NLR temperature 
range of 10^ — 2.5 x 10* K. The outward velocities obtained in 
our simulations (few 100 km/s) are smaller than the upper 
limit of few 1000 km/s observed in X-ray spectroscopic stud- 
ies of AGN warm absorber outflows (iTurner et al. 20051: 
Evans et all l2007l: [Markowitz et al.l l2008l: iTurner et all 



20081: iMehdipour. Branduardi-Ravmont fc Pag^ I2OI0I : 



Turner et al.ll2010l ). which could be due to the absence of 



radiation pressure in our current simulations. 

Thermal conduction becomes important in order to de- 
scribe the physical evol ution of TI and to resolve t he critical 
length scale of TI (e.g.. iKovama fc Inutsu"kall2004l ). The rel- 
evant scale is the Field length (Af), below which thermal 
conduction suppresses the growth of density perturbations 
due to TI by efficiently conducting heat. We compute Af for 
our cooling function, finding that it increases monotonically 
with increasing radius and Lx ■ The highest values are in the 
run with Lx ~ 0.05, which lies between Af ~ 4 x 10~® pc 
near nn, and Af ~ O.lpc near rout. The smoothing lengths 
in our simulations are ~ 0.07 pc near nn, and ~ 10 pc near 
rout. Thus our spatial resolution limit is at least 10^ — 10* 
times larger than Af. Since we do not resolve the scales 
where the growth of TI- induced density perturbations is sup- 
pressed by heat conduction, our simulations can only give 
an upper limit on the size of the smallest clumps. 

In order to check if our simulations can follow the 
isobaric mode of TI (i.e., perturbations with a constant- 
pressure), we compare the sound crossing time [ts) with 
the cooling time {tc), both local values calculated for each 
gas particle. The GADGET code determines the adaptive 
smoothing length (hsmi) of each particle to include a con- 
stant mass for the estimated density, amounting to 32 — 64 
smoothing neighbors within hsmi ■ In our simulations we used 
50 ± 1 as the desired number of SPH smoothing neighbors. 
Therefore within Shsrni there are ~150 neighboring parti- 
cles, which is a typical number needed to resolve a sound 
wave. We determine ts by computing how long it would 
take a sound wave at a given radius to cross 3/i3,„i . We find 
that ts > tc everywhere before the gas undergoes Tl-driven 
cooling, with [ts/tc]niin ~ 2. StcUirigwcrf (1982) studied the 
local and global stability of spherically symmetric. X-ray 
heated and optically thin accretion flow. He gave a criti- 
cal wavelength Ao = 4:-K'y^^^Csu/{du/dt)cooi, below which 
the instability growth time is independent of A. Here Ca 
is the sound speed, u is the specific internal energy, and 



{du/dt)cooi is the cooling term in the specific internal energy 
equation. Wavelengths shorter than Ao are most unstable, 
and such perturbations are approximately isobaric. Compar- 
ing Ao with hsmi in our simulations, we see that hsmi < Ao 
before Tl-induced cooling starts to happen. Hence we can 
resolve isobaric modes of TI. 

We check the Jeans instability criterion in order to de- 
termine if the cool, dense clumps formed in our simulations 
would undergo self-gravitating collapse. Our results indicate 
that self-gravity is sub-dominant over gas pressure. Taking 
the density and temperature of two typical clumps, we calcu- 
late the corresponding Jeans length (Aj) and mass (Mj). For 
a filament-like clump with p = 10"^° g/cm^, T = 3.2 x 10^ K, 
and of length ~ 30 pc, from Figure [T] Aj = 107.9 pc, 
and Mj = 9.7 x IO'^Mq. For a more-spherical clump with 
p = 5.0 X 10"^° g/cm^ , T = 3.2 X 10® K, and of length ~ 1 pc, 
from Figure [1 \j = 152.4 pc, and Mj = 1.4 x 10^ Mq. As 
a comparison the initial mass in the whole computational 
volume is Mtotjc = 9.77 x W^Mq (Table [T]). Therefore the 
computed Jeans length and mass are greater than the actual 
length and mass of the clumps, which means they would not 
self-gravitate, and gas pressure would restore balance over 
self-gravity. 

We compare our results to those from the Eulerian code 
ZEUS, where spherical-symmetry can be maintained to a 
higher degree. Results from ID and 2D ZEUS counterpart 
simulations show the formation of cold, dense, spherically- 
symmetric shells due to TI. These shells form continuously, 
at non-uniform time intervals, accrete in and cause large 
spikes in the central mass accretion rate. However the fiow 
remains spherical in ZEUS. We deduce that small non- 
spherical density perturbations present in SPH (from the 
initial condition), which are random, grow via TI resulting 
in the inhomogeneous fiow presented in iJ21 This has been 
tested by adding non-spherical density fiuctuations to the 
initial condition of ZEUS simulations; then the flow becomes 
clumpy and inhomogeneous in ZEUS as well. Therefore we 
cannot study the clumps in a very controlled manner using 
SPH. 

Examining the parameter space of the results described 
here, we see the occurrence of Tl-induced clumping for a rel- 
atively narrow range of X-ray luminosity (Lx = 0.01 — 0.05). 
But this is only for a fixed density (poo ~ 10^^'^ g/cm'^). At 
luminosities below Lx ~ 0.01 stable spherically-symmetric 
gas inflow occurs, and above Lx = 0.05 there is a spherical 
outflow. The range of photoionization parameter (^ = -^) 
for the operation of TI depends on Lx , r and n. We also ob- 
serve non-spherical cooling in a run with p^o = 10~^* g/cm^, 
and Lx = 5 x 10~*. This supports our assertion that such 
fragmentation corresponds to the accreting gas attaining the 
^-range for TI (largest slope of [Trad—?] curve, or dT/d'E < 0, 



KroUkI (f 19881 ) invoked a finite-amplitude TI for the ori- 



gin of broad emission line clouds in AGN. His Figure 1 shows 
the operation of TI in the [T — H] plane. However in our sim- 
ulations the initial conditions of the gas are different, there- 
fore we see a different time evolution. We rather find a typ- 
ical r — H transition of the kind shown in the bottom-right 
panel of Figure |4] and described in i|3.1.1l for the particles 
that experience TI. 

Our results show that fragmentation can occur in a sim- 
ple system with minimal physical processes: a spherically ac- 
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creting gas can cool non-spherically and become clumpy due 
to thermal instability and radiative heating/cooling. There 
is of course the additional contribution of small density fluc- 
tuations inherent to SPH. At the same time, our analyses 
using the mesh code ZEUS indicate that when a steady solu- 
tion is very unstable, it is hard to avoid the development of 
TI. The main difference with various numerical techniques is 
the amplitude and type of perturbation. Operating in such 
a manner TI can induce, or at least enhance, star-formation 
near AGN, by causing the gas to cool and fragment. Global 
TI can be a triggering mechanism for the variability ob- 
served in the luminosity of compact objects, by making the 
accreting gas inhomogeneous. 

In our simulations the gas undergoes spherical accre- 
tion initially, and the radiation from the central SMBH is 
modeled in the form of a spherical X-ray corona. Possible fu- 
ture work would include breaking such spherical symmetry, 
e.g., considering radiation from an optically thic k, geomet- 
rically thin, standard accretio n disk (jProgal 20071). including 
rotat ion of the accreting gas (|Proga. Ostr ikcr fc Kurosawa 
20081) ■ and radiation force from the disk ({Kurosawa fc Proga 
2008li. 



5 SUMMARY AND CONCLUSION 

Gas accretion onto a central SMBH was studied in Paper I 
by performing simulations of spherical Bondi accretion using 
the 3D SPH code GADGET-3. The simulations included ra- 
diative heating by a central X-ray corona and radiative cool- 
ing. All the results presented in Paper I were of spherically- 
symmetric gas properties. However, for a certain parameter 
range of X-ray luminosity Lx/L^dd = 0.01,0.02,0.05 with 
a fixed density at the outer boundary poo = 10~^'^g/cm'^ 
(i.e., within the unstable ^-range, tj3.1.1|) . we observe non- 
spherical, multi-phase gas motion. In the current work we 
present detailed analyses of the inhomogeneous gas behav- 
ior. The main results are summarized below: 

(1) The gas goes through various modes of accretion 
for different Lx- There is a stable, spherically-symmetric 
inflow with Lx = 0.005 (Run 25). The gas is heated up 
and expands as Lx is increased, which decreases the central 
mass inflow at a varying rate. For some intermediate Lx the 
accretion flow undergoes unstable motion developing non- 
spherical features. With high enough Lx a strong outflow is 
produced, ejecting most of the gas outside. 

(2) With Lx = 0.01 and 0.02, the gas cools and be- 
comes denser in a non-spherical manner, fragments, and 
forms multiple clumps. The clumping occurs over r < 30pc 
for Lx = 0.01, and r ~ 40 - 100 pc for Lx = 0.02. It starts 
to develop from ~ 1.6 — 1.7 Myr, and continues up to the 
simulation end. The clumps are filamentary in structure, su- 
personic, and have a lower photoionization. Hotter, lighter, 
subsonic gas exist in-between the filaments, and start to flow 
outward at later times. The clumps undergo a dynamic evo- 
lution: become stretched, fragment, merge with each other, 
finally accreting into rin. 

The gas particle properties versus radius show large 
scatters (few 100s to 1000) at the radial ranges where frag- 
mentation occurs, corresponding to the multi-phase struc- 
ture of the gas. The mass infiow rate shows large fluctua- 
tions in the run with Lx = 0.01 after 1.6 Myr, and in the run 



with Lx = 0.02 after 2.3 Myr, because of the accretion of 
multi-phase gas, with the spikes corresponding to the infall 
of cold dense clumps. 

(3) As they move inward the cooled clumps remain 
denser, however heat up to above 10^ K. The central heat- 
ing at r < 1 pc is mostly by adiabatic compression reach- 
ing T > lO'^ K at Tin, while radiative heating would dom- 
inate at a few pc for some particles. The photoionization 
parameter increases while moving to smaller radii, and the 
clumps remain supersonic. Artiflcial viscosity heating in- 
creases in magnitude toward rin, however it always remains 
sub-dominant compared to adiabatic heating. This is in con- 
trast to Paper I due to the different conditions here because 
of the higher Lx ■ 

(4) The cooling and clumping is caused by the inter- 
play of thermal instability in the inflowing gas acted upon 
by radiative processes, as well as being adiabatically com- 
pressed and heated. Gas particles in our simulations undergo 
TI within the unstable ^-range, where the Trad — C equilib- 
rium curve has the steepest slope (between ^ ~ 100 — 500), 
and most physically relevant processes (those with constant- 
pressure or free falling) can give rise to perturbations with 
a shallower slope. 

We detect a distinct mode of cooling/heating in the 
[T — ^] plane. A particle while being dominated by adia- 
batic process moves up to higher-T along the unstable part 
of Trad — This happens because the adiabatic free-fall re- 
lation of T versus ^ {Tff^a oc f//,a) has the same slope as 
the radiative equilibrium curve in the unstable regions. At a 
later time the particle has a transition to being dominated by 
radiative process. Then classical TI operates, and it under- 
goes catastrophic cooling or heating (whichever dominates), 
and moves away from the unstable T^ad — C part. It finally 
settles on one of the stable branches where Trad — C has a 
shallower slope. 

(5) Typically the amplitude of temperature and den- 
sity perturbations versus time shows that the perturbations 
have an exponential growth at first, and then come to a 
saturation for the remaining simulation run. The minimum 
temperature has an exponential decrease initially, indicating 
the growth of cooling instability. It remains almost constant 
after that at the temperature of the lower stable Trad — £, 
branch, where T is a weak function of f . In terms of density 
the maximum value rises indicating the growth of cooled, 
dense clumps caused by TI. 

(6) We have multi-phase medium (cold gas with 2 x 
10'* K and hot gas with 10^ — lO'^ K) at the same radius, form- 
ing self-consistently in the simulation volume. Furthermore 
we resolve the Bondi radius of the accretion flow. This can be 
used to constrain the a-parameter that is multiplied to the 
Bondi-Hoyle-Lyttleton mass accretion rate in the sub-grid 
model of SMBH feedback in large-scale galaxy simulations. 

Considering a limiting temperature of 10^ K between 
cold and hot phases, we compute the ratio (x) of mass in- 
flow rates of cold gas over hot. We find that x is a dynamical 
quantity depending on several factors: time, spatial location, 
and X-ray luminosity. The maximum value of x varies be- 
tween 1 — 5. At a fixed radius, x increases with time initially, 
and later remains almost-constant at a value x ~ 3 — 4. 
Thus we do not find a factor as large as a ~ 100, since our 
X 1 — 5. In other words, in our simulations we do not find 
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such a large fraction of ga.s hidden in the cold phase that is 
unresolved in cosmological simulations. 

(7) With Lx = 0.02 (Run 27) we see both fragmented 
cold clumps and buoyantly rising hot bubbles. Initially there 
is a central spherical outflow at r < 40 pc, which decreases 
the mass inflow rate at rin by ~ 1000 times between 1.4 — 
2.2 Myr. The surrounding inflowing gas between [40 — 70] pc 
cools and fragments by TI. The dense clumps grow with 
time, and continue to move inward, disrupting the spherical 
symmetry of the central hotter outflow. 

The clumps start accreting into nn att^ 2.2 Myr, caus- 
ing the central mass inflow rate to rise again. The gas motion 
becomes inhomogeneous with cooler, denser clumps falling 
in to the center but heating up as they move in; and hot- 
ter, lighter gas moving out. There are also few hot bubbles 
buoyantly rising from the center. 

(8) With a high enough value of the X-ray luminosity 
{Lx ~ 0.05 in Run 28), the gas is heated up strongly, cre- 
ating a spherically symmetric outflow. This expels a large 
fraction of gas out of rout, reduces the central mass inflow 
rate drastically, and produces a gas density profile that de- 
creases towards center. There is also the formation of an 
anisotropic, hot bubble buoyantly rising from the center. It 
has a narrow elongated structure, while its head has a bow 
shock-like morphology, where the interior gas collides with 
the relatively cold ambient medium and produces a shock 
region larger than the jet. The bubble has a greater out- 
ward velocity, lower density, higher temperature and higher 
entropy compared to the gas at the same radius. 
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Table 1. Simulations and Flow Features 



Run 
No. ^ 


Lx 


Spherical ^ 


Clumping ° 


Bubble ^ 


Motion 
(r < 40 pc, 1.4 - 2.3 Myr) 


ZO 


u.uuo 


les 


INO 




Inflow 


26 


0.01 


No 


Yes 


No 


Inflow 


27 


0.02 


No 


Yes 


Yes 


Outflow 


28 


0.05 


No 


No 


Yes 


Outflow 



The run numbers are from Table 2 of Paper I. The common parameter values are 
quoted below. These runs are started using the particle states in Run 23 (Paper I) at 
t = 1.4 Myr as the initial condition: pinit = pRun23, Vinit = ■yRun23, Unit = TRun23, and 
changing Lx in each case. All the runs have 7run = 5/3, nn = 0.1 pc, rout = 200 pc, 
N = 1.24 X 10^ A/totjc = 9.77 x 10'^ M©, and Mp^rt = O.791M0. The IC for Run 23 
was generated using 7init = 5/3, Poo = 10~^^ g/cm^. Too = 10^ K (for which the Bondi 
radius is Rb = 183.9 pc), pmit(r') = ps(r), ?;init(r') = vs(r), and Tinit = Trad. 

^ Whether the gas properties in the flow remain spherically-symmetric. 

Whether the gas undergoes non-uniform cooling, fragmentation and form clumps. 

Whether there is the formation of hot bubble(s) rising buoyantly from the center. 
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Figure 1. Cross-sections in Run 26 (ijf /^Edd = 0.01) showing the inner 30 pc of the [x — y] plane through z = at a time t = 2.047 Myr. 
The gas density is in the top-left panel, temperature in the top-right, photoionization parameter in the bottom-left, and Mach number 
in the bottom-right, all in cgs units. The velocity vectors are overplotted as arrows, whose sizes are denoted at the top-left of each panel 
in km/s. It shows colder, denser filament-like structures, with hotter, less-dense gas in-between, both components accreting in (with the 
colder phase moving in faster). These features are caused by Tl-induced non-spherical cooling and fragmentation. Note that a different 
color scheme has been used to depict each quantity, therefore the map between the colors and the order of values (increasing/decreasing) 
could be different from one panel to the other. This and all the other cross-section images in this paper have been generated using 
SPLASH (,Price„2007i) . 
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Figure 2. Same as Figure [T] but showing zoom-in of the inner 4pc of the [x — y] plane. It shows stretching of the colder clumps as 
they fall in toward the center. They remain denser, however heats up at r < 1 pc, mostly by adiabatic compression. Note that the color 
scheme in this cross-section has been changed and it has been plotted without the velocity vectors, in order to show the small-scale 
features clearly. 
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Figure 3. Properties of particles in Run 26 at i = 2.0 Myr (0.047 Myr earlier than in Figures □ andEll. The top two rows show four 
quantities vs. radius: radial velocity, density, temperature and Mach number. The bottom row plots the temperature vs. photo- and 
pressure-ionization parameters. Note that lower photoionization parameter (g = ^j^) corresponds to large radii, therefore the points on 
the left-hand-side in the top two rows are on the right-hand-side in the bottom-left panel. The red curve shows the free-fall scaling of 
the corresponding quantity in the top four panels, and the radiative equilibrium scaling in the bottom two panels. The blue horizontal 
straight line in the top-left and middle-right panels indicate Vr = and Mach = 1. 
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Figure 4. Time evolution of a single particle in Run 26 as it moves inward. Letters A, B, C, D, E, F, G denote the major transition 
points on the track. The starting point (A) is denoted by the triangle in each panel, showing the initial condition at i = 1.4 Myr and 
r = 53 pc, The ending point (G) at t = 1.8 Myr and r = 0.99 pc is represented by the square in each panel. The plus signs denote 
the relevant quantity in uniform time intervals of 0.004 Myr, except the last two points (inner-most in r: F, G) which are separated by 
~ 0.001 Myr. The top two rows show four radial properties: radial and angular velocities, density, temperature and entropy. The bottom 
row plots the evolution in the temperature vs. photo- and pressure-ionization parameter planes, where the direction of progress of time 
is indicated by arrows. The red curve denotes the free-fall scaling of the corresponding quantity in the top two and middle-left panels, 
and the radiative equilibrium scaling in the bottom tw o panels. The green curve in the bottom-left panel is the T — ( solution of the 
instability criterion from Equation (5) of lBalbi3 l|l98(J l. The slopes of three characteristic processes are shown as the blue lines in the 
bottom row: adiabatic free-fall (T ~ and T ~ as the dotted line, constant-pressure (T ~ ^, and S = constant) as the dashed 

line, and constant-density (^ = constant, and T ~ as the dash-dotted line. These slopes are used to describe various relevant 

perturbations (see t|3.1.1l for more details). 
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Figure 5. Ratio of thermal instability growth time (tgi) and free-fall time (iff) in the top panel at three simulation epochs: t = 1.42 Myr 
(dashed curve), 1.5 Myr (dotted), and 1.56 Myr (solid), when the flow is still spherical. The ratio is plotted only at the radial ranges with 
a positive growth rate and represent regions where TI can grow. The radial gaps indicate regions which have negative growth rate and 
are stable to TI. Mach number in the bottom panel; the upper and lower curve of each linestyle denote the median and minimum Mach 
at a given radius. Details are described in i|3.1.2l 
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Figure 6. Growth of temperature (top panel) and density (bottom panel) perturbation caused by TI in Run 26, at radius r = lOpc 
within a width Ar = Ipc, as a function of time. Minimum (dashed), maximum (solid), and average (dotted) values within the radial bin 
are plotted as a ratio of the stable unperturbed average value at t = 1.42 Myr. The relevant perturbations represented by Tmim pmax, 
and (p) show an exponential growth (linear on the log scale of the y-axis) for some time initially, and then come to a saturation. Details 
are described in i|3.1.2l last two paragraphs. 
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Figure 7. Spatial and temporal evolution of the ratio x (defined as the ratio of mass inflow rates of cold phase gas over hot, below and 
above the limiting temperature of 10^ K) in Run 26. Left panel shows the radial variation of x four simulation times: t = 1.6 Myr 
(dash-dotted curve), 1.7Myr (dashed), 1.8 Myr (solid), and 2.0 Myr (dotted), x is non-zero only within r ~ 1 — 30pc where there is cold 
gas. Right panel shows the time evolution of x at a radius r = 10 pc. At a fixed-r, x increases with time at first, and then becomes 
saturated at a value x ~ 3 — 4. Details are in ^3.1.31 
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Figure 8. Time evolution of gas in Run 27 (ijf /^Edd = 0.02) showing the inner 100 pc of the [x — 2] plane through y = 0. The gas 
density is in the left panel of each row and temperature in the right panel, all in cgs units. The velocity vectors are overplotted as arrows, 
whose sizes are denoted at the top-left of each panel in km/s. Rows correspond to the times: t = 1.86 Myr (top), t = 2.12 Myr (middle), 
and t = 2.46 Myr (bottom). A central spherically-symmetric outflow inside r ~ 40pc can be seen in the top row, with surrounding gas 
starting to cool and clump. The spherical-symmetry of the outflow is lost in the middle row, and the formation of a hot, less-dense bubble 
can be seen toward the right just off-center. The surrounding cool clumps are well-formed and fragmented, and are moving inward. The 
bottom row shows an inhomogeneous mixture of multi-phase gas motion. The denser clumps continue falling into the center, but heating 
up as they move in. The hotter gas moves outward, with signatures of few hot bubbles buoyantly rising from the center. 
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Figure 9. Time evolution of gas in Run 28 (Lx/L^dd = 0.05) showing the whole computational volume 200 pc of the [y — z] plane 
through X = 0. The gas density is in the left panel of each row and temperature in the right panel, overplotted with the velocity vector 
arrows. Rows correspond to the times: t = 1.8 Myr (top), and t = 3.0 Myr (bottom). The gas is outflowing in almost all over the volume, 
except near the very center toward the right in the plotted plane. A hot, less-dense bubble is well-formed and buoyantly rises from the 
center along the negative z-axis. The rest of the outflowing gas remains spherically symmetric. 
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Figure 10. Properties of gas in Run 28 at t = 2.0 Myr shown as the cyan and black plus symbols. The cyan points are particles 
overlapping with the negative-y axis (through their smoothing lengths) at r < 60 pc hence belonging to the bubble, and the black points 
are particles overlapping with the positive-z axis representing rest of the spherically-symmetric volume. The top two rows show four 
quantities vs. radius: radial velocity, density, temperature and entropy. The bottom row plots the temperature vs. photo- and pressure- 
ionization parameters. The red curve denotes the free-fall scaling of the corresponding quantity in the top two panels, and the radiative 
equilibrium scaling i n the bottom t wo panels. The green curve in the bottom-left panel is the T — % solution of the instability criterion 
from Equation (5) of lBalbusI l ll986l) . The blue curve overplotted in each panel represent the time evolution of a single particle inside the 
bubble, as it moves first inward and then outward, in a format similar to Figure |4] Letters A, B, C, D, E, F, G (written in blue) mark 
the major transition points of the bubble particle along the track. The starting point A from the initial condition at t = 1.4 Myr when it 
is at r = 11 pc is indicated by the triangle, the ending point G at t = 3.8 Myr when it is at r = 163 pc by the square, and the direction 
of time evolution is shown by an arrow in each panel. The orange filled circle denotes the particle at the same time t = 2.0 Myr of the 
snapshot. 



© 0000 RAS, MNRAS 000, 000-000 



24 P. Barai et al. 



> 

CD 



100.0 - 



10.0 r 



1.0 



0.1 , 



^ I I I ri\ I I I I I I If} 



1.42 (Myr)" 
1.46 - 
1.52 
1.6 



rMm 




0.1 



1.0 



10.0 

(PC) 



100.0 



Figure 11. Ratio of the Brunt- Vaisala timescale {tgy) for convective instability and free-fall time (tff) at four simulation epochs when 
the flow is nearly steady: t = 1.42Myr (dotted), 1.46Myr (dashed), 1.52Myr (solid), and 1.6Myr (dot-dash). The ratio is plotted only 
at the radial ranges having a real effective Brunt- Vaisala frequency {ij^gy > 0) and represent convectively unstable regions. The radial 



gaps indicate regions which have imaginary frequency {^%y < 0) and are convectively stable. The details are described in 
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Figure 12. Mass inflow rate A/in ^i,, across the inner boundary as a function of time for the four runs in Table[T] with Lx/^Edd = 0.005 
(Run 25), 0.01 (Run 26), 0.02 (Run 27), and 0.05 (Run 28). The mass inflow rate decreases as Lx is increased, which is associated with 
the development of an outflow. The noisy Min^.^ in Runs 26 and 27 is because of the accretion of multiphase medium, with the spikes 
corresponding to the accretion of cold dense clumps. The dash-dotted horizontal line marks the Bondi accretion rate corresponding to 
the Poo and Too of the simulations. The details are discussed in i) l3.4l 
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Figure 13. Mass outflow rate as a function of radius and time for tiie tliree runs, witii Lx = 0.01 (Run 26, top), 0.02 (Run 27, middle), 
and 0.05 (Run 28, bottom). Afout increases as Lx is increased, which is associated with the development of strong outflows in runs with 
higher Lx- The details are discussed in ii l3.4l 
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